Usenet 1993 July/InfoMagic USENET CD-ROM July 1993.ISO/sources/misc/volume15/ggems/part03

Text File | 1990-10-14 | 53.5 KB | 1,997 lines

Newsgroups: comp.sources.misc
X-UNIX-From: [e-mail]
subject: v15i025: Graphics Gems, Part 3/5
from: Craig Kolb <[e-mail]>
Sender: [e-mail] (Brandon S. Allbery - comp.sources.misc)

Posting-number: Volume 15, Issue 25
Submitted-by: Craig Kolb <[e-mail]>
Archive-name: ggems/part03

#! /bin/sh
# This is a shell archive.  Remove anything before this line, then unpack
# it by saving it into a file and typing "sh file".  To overwrite existing
# files, type "sh file -c".  You can also feed this as standard input via
# unshar, or by typing "sh <file", e.g..  If this archive is complete, you
# will see the following message at the end:
#        "End of archive 3 (of 5)."
# Contents:  AALines/utah.c Albers.c BoundSphere.c Dissolve.c
#   DoubleLine.c GraphicsGems.h LineEdge.c MatrixInvert.c
#   PolyScan/poly_clip.c PolyScan/poly_scan.c Roots3And4.c
# Wrapped by craig@weedeater on Fri Oct 12 15:53:13 1990
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f 'AALines/utah.c' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'AALines/utah.c'\"
else
echo shar: Extracting \"'AALines/utah.c'\" \(4532 characters\)
sed "s/^X//" >'AALines/utah.c' <<'END_OF_FILE'
X/*
X    file:        utah.c
X    description:    interface to Utah RLE toolkit
X    author:        A. T. Campbell
X    date:        October 27, 1989
X*/
X
X#ifndef lint
Xstatic char    sccsid[] = "%W% %G%";        /* SCCS info */
X#endif lint
X
X#include <math.h>
X#include <stdio.h>
X#ifdef sequent
X#include <strings.h>
X#else
X#include <string.h>
X#endif
X#include "utah.h"
X
X/******************************************************************************/
X
X/* return values */
Xextern void    free();
Xextern char    *malloc();
X
X/******************************************************************************/
X
Xutah_read_close(ufp)
XUTAH_FILE    *ufp;
X{
X    return(0);
X}
X
X/******************************************************************************/
X
XUTAH_FILE *
Xutah_read_init(fname, ht, wd)
X
Xchar    *fname;
Xint    *ht, *wd;
X{
X    FILE        *fp;
X    UTAH_FILE    *ufp;
X
X    /* open output stream */
X    if (!strcmp(fname, ""))
X        fp = stdin;
X    else {
X        if ((fp = fopen(fname, "r")) == NULL)
X        return(NULL);
X    }
X
X     /* change the default sv_globals struct to match what we need */
X    ufp = (UTAH_FILE *) malloc(sizeof(UTAH_FILE));
X    *ufp = sv_globals;
X    ufp->svfb_fd = fp;
X
X    /* read the header in the input file */
X      if (rle_get_setup(ufp) != 0)
X        return(NULL);
X
X    /* get image size */
X    *wd = ufp->sv_xmax - ufp->sv_xmin + 1;
X    *ht = ufp->sv_ymax - ufp->sv_ymin + 1;
X
X    /* normal termination */
X    return(ufp);
X}
X
X/******************************************************************************/
X
Xutah_read_pixels(ufp, pixels)
X
XUTAH_FILE     *ufp;
Xunsigned char    pixels[][3];
X{
X    static unsigned    n = 0;
X    static unsigned char    *r = NULL, *g = NULL, *b = NULL;
X    int        i, width;
X
X    /* allocate storage */
X    width = ufp->sv_xmax + 1;
X    if (width > n) {
X        if (n > 0) {
X            free((char *)r);
X            free((char *)g);
X            free((char *)b);
X        }
X        n = width;
X        r = (unsigned char *) malloc(n * sizeof(unsigned char));
X        g = (unsigned char *) malloc(n * sizeof(unsigned char));
X        b = (unsigned char *) malloc(n * sizeof(unsigned char));
X    }
X
X    /* read this row */
X    utah_read_rgb(ufp, r, g, b);
X
X    /* convert to pixels */
X    for (i = 0; i < width; i++) {
X        pixels[i][0] = r[i];
X        pixels[i][1] = g[i];
X        pixels[i][2] = b[i];
X    }
X
X    return(0);
X}
X
X/******************************************************************************/
X
Xutah_read_rgb(ufp, r, g, b)
X
XUTAH_FILE    *ufp;
Xunsigned char    r[], g[], b[];
X{
X    rle_pixel    *rows[3];
X
X    /* set color channels */
X    rows[0] = r;
X    rows[1] = g;
X    rows[2] = b;
X
X    /* read this row */
X    rle_getrow(ufp, rows);
X    return(0);
X}
X
X/******************************************************************************/
X
Xutah_write_close(ufp)
X
XUTAH_FILE    *ufp;
X{
X    if (!ufp) return(-1);
X    sv_puteof(ufp);
X    return(0);
X}
X
X/******************************************************************************/
X
XUTAH_FILE *
Xutah_write_init(fname, ht, wd)
X
Xchar    *fname;
Xint    ht, wd;
X{
X    FILE        *fp;
X    UTAH_FILE    *ufp;
X
X    /* open output stream */
X    if (!strcmp(fname, ""))
X        fp = stdout;
X    else {
X        if ((fp = fopen(fname, "w")) == NULL)
X        return(NULL);
X    }
X
X     /* change the default sv_globals struct to match what we need */
X    ufp = (UTAH_FILE *) malloc(sizeof(UTAH_FILE));
X    *ufp = sv_globals;
X    ufp->svfb_fd = fp;
X    ufp->sv_xmax = wd - 1;
X    ufp->sv_ymax = ht - 1;
X    ufp->sv_alpha = 0;    /* No coverage (alpha) */
X
X    /* create the header in the output file */
X      sv_setup(RUN_DISPATCH, ufp);
X
X    /* normal termination */
X    return(ufp);
X}
X
X/******************************************************************************/
X
Xutah_write_pixels(ufp, pixels)
X
XUTAH_FILE    *ufp;
Xunsigned char    pixels[][3];
X{
X    static unsigned    n = 0;
X    static unsigned char    *r = NULL, *g = NULL, *b = NULL;
X    int        i, width;
X
X    /* allocate storage */
X    width = ufp->sv_xmax + 1;
X    if (width > n) {
X        if (n > 0) {
X            free((char *)r);
X            free((char *)g);
X            free((char *)b);
X        }
X        n = width;
X        r = (unsigned char *) malloc(n * sizeof(unsigned char));
X        g = (unsigned char *) malloc(n * sizeof(unsigned char));
X        b = (unsigned char *) malloc(n * sizeof(unsigned char));
X    }
X
X    /* convert to color channels */
X    for (i = 0; i < width; i++) {
X        r[i] = pixels[i][0];
X        g[i] = pixels[i][1];
X        b[i] = pixels[i][2];
X    }
X
X    /* write this row */
X    utah_write_rgb(ufp, r, g, b);
X    return(0);
X}
X
X/******************************************************************************/
X
Xutah_write_rgb(ufp, r, g, b)
X
XUTAH_FILE    *ufp;
Xunsigned char    r[], g[], b[];
X{
X    rle_pixel    *rows[3];
X    int        width;
X
X    /* set color channels */
X    rows[0] = r;
X    rows[1] = g;
X    rows[2] = b;
X
X    /* write this row */
X    width = ufp->sv_xmax - ufp->sv_xmin + 1;
X    sv_putrow(rows, width, ufp);
X    return(0);
X}
X
X/******************************************************************************/
END_OF_FILE
if test 4532 -ne `wc -c <'AALines/utah.c'`; then
    echo shar: \"'AALines/utah.c'\" unpacked with wrong size!
fi
# end of 'AALines/utah.c'
fi
if test -f 'Albers.c' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'Albers.c'\"
else
echo shar: Extracting \"'Albers.c'\" \(3994 characters\)
sed "s/^X//" >'Albers.c' <<'END_OF_FILE'
X/*
XAlbers Equal-Area Conic Map Projection
Xby Paul Bame
Xfrom "Graphics Gems", Academic Press, 1990
X*/
X
X/*
X * Albers Conic Equal-Area Projection
X * Formulae taken from "Map Projections Used by the U.S.
X * Geological Survey" Bulletin #1532
X *
X * Equation reference numbers and some variable names taken
X * from the reference.
X * To use, call albers setup() once and then albers_project()
X * for each coordinate pair of interest.
X*/
X
X#include "GraphicsGems.h"
X#include <stdio.h>
X#include <math.h>
X
X/*
X * This is the Clarke 1866 Earth spheroid data which is often
X * used by the USGS - there are other spheroids however - see the
X * book.
X */
X
X/*
X * Earth radii in different units */
X#define CONST_EradiusKM (6378.2064)     /* Kilometers */
X#define CONST_EradiusMI (CONST_EradiusKM/1.609)    /* Miles */
X#define CONST_Ec        (0.082271854)    /* Eccentricity */
X#define CONST_Ecsq    (0.006768658)    /* Eccentricity squared */
X
X
X
X/*
X * To keep things simple, assume Earth radius is 1.0. Projected
X * coordinates (X and Y obtained from albers project ()) are
X *  dimensionless and may be multiplied by any desired radius
X *  to convert to desired units (usually Kilometers or Miles).
X*/
X#define CONST_Eradius 1.0
X
X/* Pre-computed variables */
Xstatic double middlelon;        /* longitude at center of map */
Xstatic double bigC, cone_const, r0;    /* See the reference */
X
Xstatic
Xcalc_q_msq(lat, qp, msqp)
Xdouble lat;
Xdouble *qp;
Xdouble *msqp;
X/*
X * Given latitude, calculate 'q' [eq 3-12]
X * if msqp is != NULL, m^2  [eq. 12-15].
X*/
X{
X    double s, c, es;
X
X    s = sin(lat);
X    es = s * CONST_Ec;
X
X    *qp = (1.0 - CONST_Ecsq) * ((s / (1 - es * es))-
X        (1 / (2 * CONST_Ec)) * log((1 - es) / (1 + es)));
X
X    if (msqp != NULL)
X    {
X        c = cos(lat);
X        *msqp = c * c/ (1 - es * es);
X    }
X}
X
X
X
X
Xalbers_setup(southlat, northlat, originlon, originlat)
Xdouble southlat, northlat;
Xdouble originlon;
Xdouble originlat;
X/*
X * Pre-compute a bunch of variables which are used by the
X * albers_project()
X *
X * southlat     Southern latitude for Albers projection
X * northlat    Northern latitute for Albers projection
X * originlon    Longitude for origin of projected map
X * originlat    Latitude for origin of projected map -
X *                often (northlat + southlat) / 2
X*/
X{
X    double q1, q2, q0;
X    double m1sq, m2sq;
X
X    middlelon = originlon;
X
X    cal_q_msq(southlat, &q1, &m1sq);
X    cal_q_msq(northlat, &q2, &m2sq);
X    cal_q_msq(originlat, &q0, NULL);
X
X    cone_const = (m1sq - m2sq) / (q2 - q1);
X    bigC = m1sq + cone_const * q1;
X    r0 = CONST_Eradius * sqrt(bigC - cone_const *q0) / cone_const;
X}
X
X/***************************************************************/
X
Xalbers_project(lon, lat, xp, yp)
Xdouble lon, lat;
Xdouble *xp, *yp;
X/*
X * Project lon/lat (in radians) according to albers_setup and
X * return the results via xp, yp. Units of *xp and *yp are same
X * as the units used by CONST_Eradius.
X*/
X{
X    double q, r, theta;
X    calc_q_msq(lat, &q, NULL);
X    theta = cone_const * (lon -middlelon);
X    r = CONST_Eradius * sqrt(bigC - cone_const * q) / cone_const;
X    *xp = r * sin(theta);
X    *yp = r0 - r * cos(theta);
X}
X
X#ifdef TESTPROGRAM
X
X/*
X * Test value from the USGS book. Because of roundoff, the
X * actual values are printed for visual inspection rather
X * than guessing what constitutes "close enough".
X*/
X/* Convert a degress, minutes pair to radians */
X#define DEG_MIN2RAD(degrees, minutes) \
X    ((double) ((degrees + minutes / 60.0) * M_PI / 180.0))
X
X#define Nlat DEG_MIN2RAD(29, 30)    /* 29 degrees, 30' North Lat */
X#define Slat DEG_MIN2RAD(45, 30)
X#define Originlat DEG_MIN2RAD(23, 0)
X#define Originlon DEG_MIN2RAD(-96, 0) /* West longitude is negative */
X
X#define Testlat DEG_MIN2RAD(35, 0)
X#define Testlon DEG_MIN2RAD(-75, 0)
X
X#define TestX 1885.4727
X#define TestY 1535.9250
X
Xmain()
X{
X    int i;
X    double x, y;
X
X/* Setup is also from USGS book test set */
X    albers_setup(Slat, Nlat, Originlon, Originlat);
X
X    albers_project(Testlon, Testlat, &x, &y);
X    printf("%.41f, %.41f =?= %.41f, %.41f/n",
X        x * CONST_EradiusKM, y * CONST_EradiusKM,
X        TestX, TestY);
X}
X#endif        /* TESTPROGRAM */
END_OF_FILE
if test 3994 -ne `wc -c <'Albers.c'`; then
    echo shar: \"'Albers.c'\" unpacked with wrong size!
fi
# end of 'Albers.c'
fi
if test -f 'BoundSphere.c' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'BoundSphere.c'\"
else
echo shar: Extracting \"'BoundSphere.c'\" \(3767 characters\)
sed "s/^X//" >'BoundSphere.c' <<'END_OF_FILE'
X/*
XAn Efficient Bounding Sphere
Xby Jack Ritter
Xfrom "Graphics Gems", Academic Press, 1990
X*/
X
X/* Routine to calculate tight bounding sphere over    */
X/* a set of points in 3D */
X/* This contains the routine find_bounding_sphere(), */
X/* the struct definition, and the globals used for parameters. */
X/* The abs() of all coordinates must be < BIGNUMBER */
X/* Code written by Jack Ritter and Lyle Rains. */
X
X#include <stdio.h>
X#include <math.h>
X#include "GraphicsGems.h"
X
X#define BIGNUMBER 100000000.0          /* hundred million */
X
X/* GLOBALS. These are used as input and output parameters. */
X
Xstruct Point3Struct caller_p,cen;
Xdouble rad;
Xint NUM_POINTS;
X
X/* Call with no parameters. Caller must set NUM_POINTS */
X/* before calling. */
X/* Caller must supply the routine GET_iTH_POINT(i), which loads his */
X/* ith point into the global struct caller_p. (0 <= i < NUM_POINTS). */
X/* The calling order of the points is irrelevant. */
X/* The final bounding sphere will be put into the globals */
X/* cen and rad. */
X
X
Xfind_bounding_sphere()
X{
Xregister int i;
Xregister double dx,dy,dz;
Xregister double rad_sq,xspan,yspan,zspan,maxspan;
Xdouble old_to_p,old_to_p_sq,old_to_new;
Xstruct Point3Struct xmin,xmax,ymin,ymax,zmin,zmax,dia1,dia2;
X
X
X/* FIRST PASS: find 6 minima/maxima points */
Xxmin.x=ymin.y=zmin.z= BIGNUMBER; /* initialize for min/max compare */
Xxmax.x=ymax.y=zmax.z= -BIGNUMBER;
Xfor (i=0;i<NUM_POINTS;i++)
X    {
X    GET_iTH_POINT(i); /* load global struct caller_p with */
X                    /* his ith point. */
X    if (caller_p.x<xmin.x)
X        xmin = caller_p; /* New xminimum point */
X    if (caller_p.x>xmax.x)
X        xmax = caller_p;
X    if (caller_p.y<ymin.y)
X        ymin = caller_p;
X    if (caller_p.y>ymax.y)
X        ymax = caller_p;
X    if (caller_p.z<zmin.z)
X        zmin = caller_p;
X    if (caller_p.z>zmax.z)
X        zmax = caller_p;
X    }
X/* Set xspan = distance between the 2 points xmin & xmax (squared) */
Xdx = xmax.x - xmin.x;
Xdy = xmax.y - xmin.y;
Xdz = xmax.z - xmin.z;
Xxspan = dx*dx + dy*dy + dz*dz;
X
X/* Same for y & z spans */
Xdx = ymax.x - ymin.x;
Xdy = ymax.y - ymin.y;
Xdz = ymax.z - ymin.z;
Xyspan = dx*dx + dy*dy + dz*dz;
X
Xdx = zmax.x - zmin.x;
Xdy = zmax.y - zmin.y;
Xdz = zmax.z - zmin.z;
Xzspan = dx*dx + dy*dy + dz*dz;
X
X/* Set points dia1 & dia2 to the maximally seperated pair */
Xdia1 = xmin; dia2 = xmax; /* assume xspan biggest */
Xmaxspan = xspan;
Xif (yspan>maxspan)
X    {
X    maxspan = yspan;
X    dia1 = ymin; dia2 = ymax;
X    }
Xif (zspan>maxspan)
X    {
X    dia1 = zmin; dia2 = zmax;
X    }
X
X
X/* dia1,dia2 is a diameter of initial sphere */
X/* calc initial center */
Xcen.x = (dia1.x+dia2.x)/2.0;
Xcen.y = (dia1.y+dia2.y)/2.0;
Xcen.z = (dia1.z+dia2.z)/2.0;
X/* calculate initial radius**2 and radius */
Xdx = dia2.x-cen.x; /* x componant of radius vector */
Xdy = dia2.y-cen.y; /* y componant of radius vector */
Xdz = dia2.z-cen.z; /* z componant of radius vector */
Xrad_sq = dx*dx + dy*dy + dz*dz;
Xrad = sqrt(rad_sq);
X
X/* SECOND PASS: increment current sphere */
X
Xfor (i=0;i<NUM_POINTS;i++)
X    {
X    GET_iTH_POINT(i); /* load global struct caller_p  */
X                    /* with his ith point. */
X    dx = caller_p.x-cen.x;
X    dy = caller_p.y-cen.y;
X    dz = caller_p.z-cen.z;
X    old_to_p_sq = dx*dx + dy*dy + dz*dz;
X    if (old_to_p_sq > rad_sq)     /* do r**2 test first */
X        {     /* this point is outside of current sphere */
X        old_to_p = sqrt(old_to_p_sq);
X        /* calc radius of new sphere */
X        rad = (rad + old_to_p) / 2.0;
X        rad_sq = rad*rad;     /* for next r**2 compare */
X        old_to_new = old_to_p - rad;
X        /* calc center of new sphere */
X        cen.x = (rad*cen.x + old_to_new*caller_p.x) / old_to_p;
X        cen.y = (rad*cen.y + old_to_new*caller_p.y) / old_to_p;
X        cen.z = (rad*cen.z + old_to_new*caller_p.z) / old_to_p;
X        /* Suppress if desired */
X        printf("\n New sphere: cen,rad = %f %f %f   %f",
X            cen.x,cen.y,cen.z, rad);
X        }
X    }
X}              /* end of find_bounding_sphere()  */
END_OF_FILE
if test 3767 -ne `wc -c <'BoundSphere.c'`; then
    echo shar: \"'BoundSphere.c'\" unpacked with wrong size!
fi
# end of 'BoundSphere.c'
fi
if test -f 'Dissolve.c' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'Dissolve.c'\"
else
echo shar: Extracting \"'Dissolve.c'\" \(4596 characters\)
sed "s/^X//" >'Dissolve.c' <<'END_OF_FILE'
X/* A Digital Dissolve Effect
Xby Mike Morton
Xfrom "Graphics Gems", Academic Press, 1990
X*/
X
X/*
X * Code fragment to advance from one element to the next.
X *
X * int reg;                /* current sequence element
X * reg = 1;                /* start in any non-zero state
X * if (reg & 1)                /* is the bottom bit set?
X *     reg = (reg >>1) ^ MASK;        /* yes: toss out 1 bit; XOR in mask
X * else reg = reg >>1;            /* no: toss out 0 bit
X */
X
Xint randmasks[32];    /* Gotta fill this in yourself. */
X
Xdissolve1 (height, width)    /* first version of the dissolve                                 /* algorithm */
X    int height, width;    /* number of rows, columns */
X{
X    int pixels, lastnum;    /* number of pixels; */
X                /* last pixel's number */
X    int regwidth;        /* "width" of sequence generator */
X    register long mask;    /* mask to XOR with to*/
X                    /* create sequence */
X    register unsigned long element;
X                    /* one element of random sequence */
X    register int row, column;
X                    /* row and column numbers for a pixel */
X
X      /* Find smallest register which produces enough pixel numbers */
X     pixels = height * width; /* compute number of pixels */
X                            /* to dissolve */
X     lastnum = pixels-1;    /* find last element (they go 0..lastnum) */
X     regwidth = bitwidth (lastnum); /* how wide must the */
X                    /* register be? */
X     mask = randmasks [regwidth];    /* which mask is for that width? */
X
X     /* Now cycle through all sequence elements. */
X
X      element = 1;    /* 1st element (could be any nonzero) */
X
X
X      do {
X        row = element / width;    /* how many rows down is this pixel? */
X        column = element % width;    /* and how many columns across? */
X        if (row < height)    /* is this seq element in the array? */
X          copy (row, column);    /* yes: copy the (r,c)'th pixel */
X
X        /* Compute the next sequence element */
X        if (element & 1)        /* is the low bit set? */
X          element = (element >>1)^mask;    /* yes: shift value, */
X                        /* XOR in mask */
X        else element = (element >>1);    /* no: just shift the value */
X     } while (element != 1);        /* loop until we return  */
X                        /* to original element */
X     copy (0, 0);        /* kludge: the loop doesn't produce (0,0) */
X}                        /* end of dissolve1() */
X
X
X
Xint bitwidth (N)    /* find "bit-width" needed to represent N */
X    unsigned int N;    /* number to compute the width of */
X{
X     int width = 0;    /* initially, no bits needed to represent N */
X     while (N != 0) {    /* loop 'til N has been whittled down to 0 */
X        N >>= 1;        /* shift N right 1 bit (NB: N is unsigned) */
X        width++;        /* and remember how wide N is */
X      }            /* end of loop shrinking N down to nothing *
X      return (width);    /* return bit positions counted */
X
X}                        /* end of bitwidth() */
X
X
X
Xdissolve2 (height, width)    /* fast version of the dissolve algorithm */
X    int height, width;    /* number of rows, columns */
X{
X    int rwidth, cwidth;    /* bit width for rows, for columns */
X    int regwidth;        /* "width" of sequence generator */
X    register long mask;    /* mask to XOR with to create sequence */
X    register int rowshift;    /* shift distance to get row  */
X                            /* from element */
X    register int colmask; /* mask to extract column from element */
X    register unsigned long element; /* one element of random */                                     /* sequence */
X    register int row, column;    /* row and column for one pixel */
X
X
X      /* Find the mask to produce all rows and columns. */
X
X    rwidth = bitwidth (height); /* how many bits needed for height? */
X    cwidth = bitwidth (width);  /* how many bits needed for width? */
X    regwidth = rwidth + cwidth; /* how wide must the register be? */
X    mask = randmasks [regwidth]; /* which mask is for that width? */
X
X /* Find values to extract row and col numbers from each element. */
X    rowshift = cwidth; /* find dist to shift to get top bits (row) */
X    colmask = (1<<cwidth)-1;    /* find mask to extract  */
X                        /* bottom bits (col) */
X
X      /* Now cycle through all sequence elements. */
X
X    element = 1;    /* 1st element (could be any nonzero) */
X    do {
X        row = element >> rowshift; /* find row number for this pixel */
X        column = element & colmask; /* and how many columns across? */
X        if ((row < height)    /* does element fall in the array? */
X            && (column < width)) /* ...must check row AND column */
X        copy (row, column); /* in bounds: copy the (r,c)'th pixel */
X
X        /* Compute the next sequence element */
X        if (element & 1)        /* is the low bit set? */
X        element = (element >>1)^mask; /* yes: shift value, /*
X                        /* XOR in mask */
X        else element = (element >>1); /* no: just shift the value */
X    } while (element != 1);     /* loop until we return to */
X                    /*  original element */
X
X    copy (0, 0);        /* kludge: element never comes up zero */
X}                    /* end of dissolve2() */
END_OF_FILE
if test 4596 -ne `wc -c <'Dissolve.c'`; then
    echo shar: \"'Dissolve.c'\" unpacked with wrong size!
fi
# end of 'Dissolve.c'
fi
if test -f 'DoubleLine.c' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'DoubleLine.c'\"
else
echo shar: Extracting \"'DoubleLine.c'\" \(4647 characters\)
sed "s/^X//" >'DoubleLine.c' <<'END_OF_FILE'
X/*
XSymmetric Double Step Line Algorithm
Xby Brian Wyvill
Xfrom "Graphics Gems", Academic Press, 1990
X*/
X
X#define swap(a,b)           {a^=b; b^=a; a^=b;}
X#define absolute(i,j,k)     ( (i-j)*(k = ( (i-j)<0 ? -1 : 1)))
Xint
Xsymwuline(a1, b1, a2, b2) int a1, b1, a2, b2;
X{
X    int           dx, dy, incr1, incr2, D, x, y, xend, c, pixels_left;
X    int           x1, y1;
X    int           sign_x, sign_y, step, reverse, i;
X
X    dx = absolute(a2, a1, sign_x);
X    dy = absolute(b2, b1, sign_y);
X    /* decide increment sign by the slope sign */
X    if (sign_x == sign_y)
X        step = 1;
X    else
X        step = -1;
X
X    if (dy > dx) {        /* chooses axis of greatest movement (make
X                          * dx) */
X        swap(a1, b1);
X        swap(a2, b2);
X        swap(dx, dy);
X        reverse = 1;
X    } else
X        reverse = 0;
X    /* note error check for dx==0 should be included here */
X    if (a1 > a2) {        /* start from the smaller coordinate */
X        x = a2;
X        y = b2;
X        x1 = a1;
X        y1 = b1;
X    } else {
X        x = a1;
X        y = b1;
X        x1 = a2;
X        y1 = b2;
X    }
X
X
X    /* Note dx=n implies 0 - n or (dx+1) pixels to be set */
X    /* Go round loop dx/4 times then plot last 0,1,2 or 3 pixels */
X    /* In fact (dx-1)/4 as 2 pixels are already plottted */
X    xend = (dx - 1) / 4;
X    pixels_left = (dx - 1) % 4;    /* number of pixels left over at the
X                                  * end */
X    plot(x, y, reverse);
X    plot(x1, y1, reverse);    /* plot first two points */
X    incr2 = 4 * dy - 2 * dx;
X    if (incr2 < 0) {    /* slope less than 1/2 */
X        c = 2 * dy;
X        incr1 = 2 * c;
X        D = incr1 - dx;
X
X        for (i = 0; i < xend; i++) {    /* plotting loop */
X            ++x;
X            --x1;
X            if (D < 0) {
X                              /* pattern 1 forwards */
X                plot(x, y, reverse);
X                plot(++x, y, reverse);
X                                /* pattern 1 backwards */
X                plot(x1, y1, reverse);
X                plot(--x1, y1, reverse);
X                D += incr1;
X            } else {
X                if (D < c) {
X                    /* pattern 2 forwards */
X                    plot(x, y, reverse);
X                    plot(++x, y += step, reverse);
X                    /* pattern 2 backwards */
X                    plot(x1, y1, reverse);
X                    plot(--x1, y1 -= step, reverse);
X                } else {
X                        /* pattern 3 forwards */
X                    plot(x, y += step, reverse);
X                    plot(++x, y, reverse);
X                    /* pattern 3 backwards */
X                    plot(x1, y1 -= step, reverse);
X                    plot(--x1, y1, reverse);
X                }
X                D + = incr2;
X            }
X        }        /* end for */
X
X        /* plot last pattern */
X        if (pixels_left) {
X            if (D < 0) {
X                plot(++x, y, reverse);    /* pattern 1 */
X                if (pixels_left > 1)
X                    plot(++x, y, reverse);
X                if (pixels_left > 2)
X                    plot(--x1, y1, reverse);
X            } else {
X                if (D < c) {
X                    plot(++x, y, reverse);    /* pattern 2  */
X                    if (pixels_left > 1)
X                        plot(++x, y += step, reverse);
X                    if (pixels_left > 2)
X                        plot(--x1, y1, reverse);
X                } else {
X                  /* pattern 3 */
X                    plot(++x, y += step, reverse);
X                    if (pixels_left > 1)
X                        plot(++x, y, reverse);
X                    if (pixels_left > 2)
X                        plot(--x1, y1 -= step, reverse);
X                }
X            }
X        }        /* end if pixels_left */
X    }
X    /* end slope < 1/2 */
X    else {            /* slope greater than 1/2 */
X        c = 2 * (dy - dx);
X        incr1 = 2 * c;
X        D = incr1 + dx;
X        for (i = 0; i < xend; i++) {
X            ++x;
X            --x1;
X            if (D > 0) {
X              /* pattern 4 forwards */
X                plot(x, y += step, reverse);
X                plot(++x, y += step, reverse);
X              /* pattern 4 backwards */
X                plot(x1, y1 -= step, reverse);
X                plot(--x1, y1 -= step, reverse);
X                D += incr1;
X            } else {
X                if (D < c) {
X                  /* pattern 2 forwards */
X                    plot(x, y, reverse);
X                    plot(++x, y += step, reverse);
X
X                   /* pattern 2 backwards */
X                    plot(x1, y1, reverse);
X                    plot(--x1, y1 -= step, reverse);
X                } else {
X                  /* pattern 3 forwards */
X                    plot(x, y += step, reverse);
X                    plot(++x, y, reverse);
X                  /* pattern 3 backwards */
X                    plot(x1, y1 -= step, reverse);
X                    plot(--x1, y1, reverse);
X                }
X                D += incr2;
X            }
X        }        /* end for */
X        /* plot last pattern */
X        if (pixels_left) {
X            if (D > 0) {
X                plot(++x, y += step, reverse);    /* pattern 4 */
X                if (pixels_left > 1)
X                    plot(++x, y += step, reverse);
X                if (pixels_left > 2)
X                    plot(--x1, y1 -= step, reverse);
X            } else {
X                if (D < c) {
X                    plot(++x, y, reverse);    /* pattern 2  */
X                    if (pixels_left > 1)
X                        plot(++x, y += step, reverse);
X                    if (pixels_left > 2)
X                        plot(--x1, y1, reverse);
X                } else {
X                  /* pattern 3 */
X                    plot(++x, y += step, reverse);
X                    if (pixels_left > 1)
X                        plot(++x, y, reverse);
X                    if (pixels_left > 2) {
X                        if (D > c) /* step 3 */
X                           plot(--x1, y1 -= step, reverse);
X                        else /* step 2 */
X                            plot(--x1, y1, reverse);
X                                 }
X                }
X            }
X        }
X    }
X}
X/* non-zero flag indicates the pixels needing swap back. */
Xplot(x, y, flag) int x, y, flag;
X{
X    if (flag)
X        setpixel(y, x);
X    else
X        setpixel(x, y);
X}
X
END_OF_FILE
if test 4647 -ne `wc -c <'DoubleLine.c'`; then
    echo shar: \"'DoubleLine.c'\" unpacked with wrong size!
fi
# end of 'DoubleLine.c'
fi
if test -f 'GraphicsGems.h' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'GraphicsGems.h'\"
else
echo shar: Extracting \"'GraphicsGems.h'\" \(4049 characters\)
sed "s/^X//" >'GraphicsGems.h' <<'END_OF_FILE'
X/*
X * GraphicsGems.h
X * Version 1.0 - Andrew Glassner
X * from "Graphics Gems", Academic Press, 1990
X */
X
X#ifndef GG_H
X
X#define GG_H 1
X
X/*********************/
X/* 2d geometry types */
X/*********************/
X
Xtypedef struct Point2Struct {    /* 2d point */
X    double x, y;
X    } Point2;
Xtypedef Point2 Vector2;
X
Xtypedef struct IntPoint2Struct {    /* 2d integer point */
X    int x, y;
X    } IntPoint2;
X
Xtypedef struct Matrix3Struct {    /* 3-by-3 matrix */
X    double element[3][3];
X    } Matrix3;
X
Xtypedef struct Box2dStruct {        /* 2d box */
X    Point2 min, max;
X    } Box2;
X
X
X/*********************/
X/* 3d geometry types */
X/*********************/
X
Xtypedef struct Point3Struct {    /* 3d point */
X    double x, y, z;
X    } Point3;
Xtypedef Point3 Vector3;
X
Xtypedef struct IntPoint3Struct {    /* 3d integer point */
X    int x, y, z;
X    } IntPoint3;
X
X
Xtypedef struct Matrix4Struct {    /* 4-by-4 matrix */
X    double element[4][4];
X    } Matrix4;
X
Xtypedef struct Box3dStruct {        /* 3d box */
X    Point3 min, max;
X    } Box3;
X
X
X
X/***********************/
X/* one-argument macros */
X/***********************/
X
X/* absolute value of a */
X#define ABS(a)        (((a)<0) ? -(a) : (a))
X
X/* round a to nearest integer towards 0 */
X#define FLOOR(a)        ((a)>0 ? (int)(a) : -(int)(-a))
X
X/* round a to nearest integer away from 0 */
X#define CEILING(a) \
X((a)==(int)(a) ? (a) : (a)>0 ? 1+(int)(a) : -(1+(int)(-a)))
X
X/* round a to nearest int */
X#define ROUND(a)    ((a)>0 ? (int)(a+0.5) : -(int)(0.5-a))
X
X/* take sign of a, either -1, 0, or 1 */
X#define ZSGN(a)        (((a)<0) ? -1 : (a)>0 ? 1 : 0)
X
X/* take binary sign of a, either -1, or 1 if >= 0 */
X#define SGN(a)        (((a)<0) ? -1 : 0)
X
X/* shout if something that should be true isn't */
X#define ASSERT(x) \
Xif (!(x)) fprintf(stderr," Assert failed: x\n");
X
X/* square a */
X#define SQR(a)        ((a)*(a))
X
X
X/***********************/
X/* two-argument macros */
X/***********************/
X
X/* find minimum of a and b */
X#define MIN(a,b)    (((a)<(b))?(a):(b))
X
X/* find maximum of a and b */
X#define MAX(a,b)    (((a)>(b))?(a):(b))
X
X/* swap a and b (see Gem by Wyvill) */
X#define SWAP(a,b)    { a^=b; b^=a; a^=b; }
X
X/* linear interpolation from l (when a=0) to h (when a=1)*/
X/* (equal to (a*h)+((1-a)*l) */
X#define LERP(a,l,h)    ((l)+(((h)-(l))*(a)))
X
X/* clamp the input to the specified range */
X#define CLAMP(v,l,h)    ((v)<(l) ? (l) : (v) > (h) ? (h) : v)
X
X
X/****************************/
X/* memory allocation macros */
X/****************************/
X
X/* create a new instance of a structure (see Gem by Hultquist) */
X#define NEWSTRUCT(x)    (struct x *)(malloc((unsigned)sizeof(struct x)))
X
X/* create a new instance of a type */
X#define NEWTYPE(x)    (x *)(malloc((unsigned)sizeof(x)))
X
X
X/********************/
X/* useful constants */
X/********************/
X
X#define PI        3.141592    /* the venerable pi */
X#define PITIMES2    6.283185    /* 2 * pi */
X#define PIOVER2        1.570796    /* pi / 2 */
X#define E        2.718282    /* the venerable e */
X#define SQRT2        1.414214    /* sqrt(2) */
X#define SQRT3        1.732051    /* sqrt(3) */
X#define GOLDEN        1.618034    /* the golden ratio */
X#define DTOR        0.017453    /* convert degrees to radians */
X#define RTOD        57.29578    /* convert radians to degrees */
X
X
X/************/
X/* booleans */
X/************/
X
X#define TRUE        1
X#define FALSE        0
X#define ON        1
X#define OFF         0
Xtypedef int boolean;            /* boolean data type */
Xtypedef boolean flag;            /* flag data type */
X
Xextern double V2SquaredLength(), V2Length();
Xextern double V2Dot(), V2DistanceBetween2Points();
Xextern Vector2 *V2Negate(), *V2Normalize(), *V2Scale(), *V2Add(), *V2Sub();
Xextern Vector2 *V2Lerp(), *V2Combine(), *V2Mul(), *V2MakePerpendicular();
Xextern Vector2 *V2New(), *V2Duplicate();
Xextern Point2 *V2MulPointByMatrix();
Xextern Matrix3 *V2MatMul();
X
Xextern double V3SquaredLength(), V3Length();
Xextern double V3Dot(), V3DistanceBetween2Points();
Xextern Vector3 *V3Normalize(), *V3Scale(), *V3Add(), *V3Sub();
Xextern Vector3 *V3Lerp(), *V3Combine(), *V3Mul(), *V3Cross();
Xextern Vector3 *V3New(), *V3Duplicate();
Xextern Point3 *V3MulPointByMatrix();
Xextern Matrix4 *V3MatMul();
X
Xextern double RegulaFalsi(), NewtonRaphson(), findroot();
X
X#endif
END_OF_FILE
if test 4049 -ne `wc -c <'GraphicsGems.h'`; then
    echo shar: \"'GraphicsGems.h'\" unpacked with wrong size!
fi
# end of 'GraphicsGems.h'
fi
if test -f 'LineEdge.c' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'LineEdge.c'\"
else
echo shar: Extracting \"'LineEdge.c'\" \(3852 characters\)
sed "s/^X//" >'LineEdge.c' <<'END_OF_FILE'
X/*
XFast Line-Edge Intersections on a Uniform Grid
Xby Andrew Shapira
Xfrom "Graphics Gems", Academic Press, 1990
X*/
X
X#include "GraphicsGems.h"
X
X#define OCTANT(f1, f2, f3, f4, f5, i1, s1, r1, r2) \
X    for (f1, f2, f3, nr = 0; f4; f5) { \
X        if (nr < liconst) { \
X            if (i1) \
X                r1(&C); \
X            else \
X                vertex(&C); \
X        } \
X        else { \
X            s1; \
X            if (nr -= liconst) { \
X                r2(&C); \
X                r1(&C); \
X            } \
X            else \
X                vertex(&C); \
X        } \
X    }
X
Xfind_intersections(Pptr, Qptr)
XIntPoint2   *Pptr, *Qptr;       /* P and Q as described in gem text */
X{
X    IntPoint2   P, Q;           /* P and Q, dereferenced for speed */
X    IntPoint2   C;              /* current grid point */
X    int         nr;             /* remainder */
X    int         deltax, deltay; /* Q.x - P.x, Q.y - P.y */
X    int         liconst;          /* loop-invariant constant */
X
X    P.x = Pptr->x;
X    P.y = Pptr->y;
X    Q.x = Qptr->x;
X    Q.y = Qptr->y;
X    deltax = Q.x - P.x;
X    deltay = Q.y - P.y;
X
X
X    /* for reference purposes, let theta be the angle from P to Q */
X
X    if ((deltax >= 0) && (deltay >= 0) && (deltay < deltax))
X                        /* 0 <= theta < 45 */
X        OCTANT(C.x = P.x + 1, C.y = P.y, liconst = deltax - deltay,
X            C.x < Q.x, C.x++, nr += deltay, C.y++, up, left)
X    else if ((deltax > 0) && (deltay >= 0) && (deltay >= deltax))
X                       /* 45 <= theta < 90 */
X        OCTANT(C.y = P.y + 1, C.x = P.x, liconst = deltay - deltax,
X             C.y < Q.y, C.y++, nr += deltax, C.x++, right, down)
X    else if ((deltax <= 0) && (deltay >= 0) && (deltay > -deltax))
X                       /* 90 <= theta < 135 */
X        OCTANT(C.y = P.y + 1, C.x = P.x, liconst = deltay + deltax,
X             C.y < Q.y, C.y++, nr -= deltax, C.x--, left, down)
X    else if ((deltax <= 0) && (deltay > 0) && (deltay <= -deltax))
X                       /* 135 <= theta < 180 */
X        OCTANT(C.x = P.x - 1, C.y = P.y, liconst = -deltax - deltay,
X             C.x > Q.x, C.x--, nr += deltay, C.y++, up, right)
X    else if ((deltax <= 0) && (deltay <= 0) && (deltay > deltax))
X                      /* 180 <= theta < 225 */
X        OCTANT(C.x = P.x - 1, C.y = P.y, liconst = -deltax + deltay,
X            C.x > Q.x, C.x--, nr -= deltay, C.y--, down, right)
X    else if ((deltax < 0) && (deltay <= 0) && (deltay <= deltax))
X                      /* 225 <= theta < 270 */
X        OCTANT(C.y = P.y - 1, C.x = P.x, liconst = -deltay + deltax,
X            C.y > Q.y, C.y--, nr -= deltax, C.x--, left, up)
X    else if ((deltax >= 0) && (deltay <= 0) && (-deltay > deltax))
X                     /* 270 <= theta < 315 */
X        OCTANT(C.y = P.y - 1, C.x = P.x, liconst = -deltay - deltax,
X            C.y > Q.y, C.y--, nr += deltax, C.x++, right, up)
X    else if ((deltax >= 0) && (deltay < 0) && (-deltay <= deltax))
X                     /* 315 <= theta < 360 */
X        OCTANT(C.x = P.x + 1, C.y = P.y, liconst = deltax + deltay,
X             C.x < Q.x, C.x++, nr -= deltay, C.y--, down, left)
X    else {}                    /* P = Q */
X}
X
Xvertex(I)
XIntPoint2   *I;
X{
X    /* Note: replace printf with code to process vertex, if desired */
X    (void) printf("vertex at %d %d\n", I->x, I->y);
X}
X
Xleft(I)
XIntPoint2   *I;
X{
X
X    /* Note: replace printf with code to process leftward */
X    /* intersection, if desired */
X    (void) printf("left from %d %d\n", I->x, I->y);
X}
X
Xup(I)
XIntPoint2   *I;
X{
X    /* Note: replace printf with code to process upward */
X    /* intersection, if desired */
X    (void) printf("up from %d %d\n", I->x, I->y);
X}
X
Xright(I)
XIntPoint2   *I;
X{
X    /* Note: replace printf with code to process rightward */
X    /* intersection, if desired */
X    (void) printf("right from %d %d\n", I->x, I->y);
X}
X
Xdown(I)
XIntPoint2   *I;
X{
X    /* Note: replace printf with code to process downward */
X    /* intersection, if desired */
X    (void) printf("down from %d %d\n", I->x, I->y);
X}
END_OF_FILE
if test 3852 -ne `wc -c <'LineEdge.c'`; then
    echo shar: \"'LineEdge.c'\" unpacked with wrong size!
fi
# end of 'LineEdge.c'
fi
if test -f 'MatrixInvert.c' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'MatrixInvert.c'\"
else
echo shar: Extracting \"'MatrixInvert.c'\" \(4893 characters\)
sed "s/^X//" >'MatrixInvert.c' <<'END_OF_FILE'
X/*
XMatrix Inversion
Xby Richard Carling
Xfrom "Graphics Gems", Academic Press, 1990
X*/
X
X#define SMALL_NUMBER    1.e-8
X/*
X *   inverse( original_matrix, inverse_matrix )
X *
X *    calculate the inverse of a 4x4 matrix
X *
X *     -1
X *     A  = ___1__ adjoint A
X *         det A
X */
X
X#include "GraphicsGems.h"
Xinverse( in, out ) Matrix4 *in, *out;
X{
X    int i, j;
X    double det, det4x4();
X
X    /* calculate the adjoint matrix */
X
X    adjoint( in, out );
X
X    /*  calculate the 4x4 determinent
X     *  if the determinent is zero,
X     *  then the inverse matrix is not unique.
X     */
X
X    det = det4x4( in );
X
X    if ( fabs( det ) < SMALL_NUMBER) {
X        printf("Non-singular matrix, no inverse!\n");
X        exit();
X    }
X
X    /* scale the adjoint matrix to get the inverse */
X
X    for (i=0; i<4; i++)
X        for(j=0; j<4; j++)
X        out->element[i][j] = out->element[i][j] / det;
X}
X
X
X/*
X *   adjoint( original_matrix, inverse_matrix )
X *
X *     calculate the adjoint of a 4x4 matrix
X *
X *      Let  a   denote the minor determinant of matrix A obtained by
X *           ij
X *
X *      deleting the ith row and jth column from A.
X *
X *                    i+j
X *     Let  b   = (-1)    a
X *          ij            ji
X *
X *    The matrix B = (b  ) is the adjoint of A
X *                     ij
X */
X
Xadjoint( in, out ) Matrix4 *in; Matrix4 *out;
X{
X    double a1, a2, a3, a4, b1, b2, b3, b4;
X    double c1, c2, c3, c4, d1, d2, d3, d4;
X    double det3x3();
X
X    /* assign to individual variable names to aid  */
X    /* selecting correct values  */
X
X    a1 = in->element[0][0]; b1 = in->element[0][1];
X    c1 = in->element[0][2]; d1 = in->element[0][3];
X
X    a2 = in->element[1][0]; b2 = in->element[1][1];
X    c2 = in->element[1][2]; d2 = in->element[1][3];
X
X    a3 = in->element[2][0]; b3 = in->element[2][1];
X    c3 = in->element[2][2]; d3 = in->element[2][3];
X
X    a4 = in->element[3][0]; b4 = in->element[3][1];
X    c4 = in->element[3][2]; d4 = in->element[3][3];
X
X
X    /* row column labeling reversed since we transpose rows & columns */
X
X    out->element[0][0]  =   det3x3( b2, b3, b4, c2, c3, c4, d2, d3, d4);
X    out->element[1][0]  = - det3x3( a2, a3, a4, c2, c3, c4, d2, d3, d4);
X    out->element[2][0]  =   det3x3( a2, a3, a4, b2, b3, b4, d2, d3, d4);
X    out->element[3][0]  = - det3x3( a2, a3, a4, b2, b3, b4, c2, c3, c4);
X
X    out->element[0][1]  = - det3x3( b1, b3, b4, c1, c3, c4, d1, d3, d4);
X    out->element[1][1]  =   det3x3( a1, a3, a4, c1, c3, c4, d1, d3, d4);
X    out->element[2][1]  = - det3x3( a1, a3, a4, b1, b3, b4, d1, d3, d4);
X    out->element[3][1]  =   det3x3( a1, a3, a4, b1, b3, b4, c1, c3, c4);
X
X    out->element[0][2]  =   det3x3( b1, b2, b4, c1, c2, c4, d1, d2, d4);
X    out->element[1][2]  = - det3x3( a1, a2, a4, c1, c2, c4, d1, d2, d4);
X    out->element[2][2]  =   det3x3( a1, a2, a4, b1, b2, b4, d1, d2, d4);
X    out->element[3][2]  = - det3x3( a1, a2, a4, b1, b2, b4, c1, c2, c4);
X
X    out->element[0][3]  = - det3x3( b1, b2, b3, c1, c2, c3, d1, d2, d3);
X    out->element[1][3]  =   det3x3( a1, a2, a3, c1, c2, c3, d1, d2, d3);
X    out->element[2][3]  = - det3x3( a1, a2, a3, b1, b2, b3, d1, d2, d3);
X    out->element[3][3]  =   det3x3( a1, a2, a3, b1, b2, b3, c1, c2, c3);
X}
X/*
X * double = det4x4( matrix )
X *
X * calculate the determinent of a 4x4 matrix.
X */
Xdouble det4x4( m ) Matrix4 *m;
X{
X    double ans;
X    double a1, a2, a3, a4, b1, b2, b3, b4, c1, c2, c3, c4, d1, d2, d3,             d4;
X    double det3x3();
X
X    /* assign to individual variable names to aid selecting */
X    /*  correct elements */
X
X    a1 = m->element[0][0]; b1 = m->element[0][1];
X    c1 = m->element[0][2]; d1 = m->element[0][3];
X
X    a2 = m->element[1][0]; b2 = m->element[1][1];
X    c2 = m->element[1][2]; d2 = m->element[1][3];
X
X    a3 = m->element[2][0]; b3 = m->element[2][1];
X    c3 = m->element[2][2]; d3 = m->element[2][3];
X
X    a4 = m->element[3][0]; b4 = m->element[3][1];
X    c4 = m->element[3][2]; d4 = m->element[3][3];
X
X    ans = a1 * det3x3( b2, b3, b4, c2, c3, c4, d2, d3, d4)
X        - b1 * det3x3( a2, a3, a4, c2, c3, c4, d2, d3, d4)
X        + c1 * det3x3( a2, a3, a4, b2, b3, b4, d2, d3, d4)
X        - d1 * det3x3( a2, a3, a4, b2, b3, b4, c2, c3, c4);
X    return ans;
X}
X
X
X
X/*
X * double = det3x3(  a1, a2, a3, b1, b2, b3, c1, c2, c3 )
X *
X * calculate the determinent of a 3x3 matrix
X * in the form
X *
X *     | a1,  b1,  c1 |
X *     | a2,  b2,  c2 |
X *     | a3,  b3,  c3 |
X */
X
Xdouble det3x3( a1, a2, a3, b1, b2, b3, c1, c2, c3 )
Xdouble a1, a2, a3, b1, b2, b3, c1, c2, c3;
X{
X    double ans;
X    double det2x2();
X
X    ans = a1 * det2x2( b2, b3, c2, c3 )
X        - b1 * det2x2( a2, a3, c2, c3 )
X        + c1 * det2x2( a2, a3, b2, b3 );
X    return ans;
X}
X
X/*
X * double = det2x2( double a, double b, double c, double d )
X *
X * calculate the determinent of a 2x2 matrix.
X */
X
Xdouble det2x2( a, b, c, d)
Xdouble a, b, c, d;
X{
X    double ans;
X    ans = a * d - b * c;
X    return ans;
X}
X
X
END_OF_FILE
if test 4893 -ne `wc -c <'MatrixInvert.c'`; then
    echo shar: \"'MatrixInvert.c'\" unpacked with wrong size!
fi
# end of 'MatrixInvert.c'
fi
if test -f 'PolyScan/poly_clip.c' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'PolyScan/poly_clip.c'\"
else
echo shar: Extracting \"'PolyScan/poly_clip.c'\" \(4217 characters\)
sed "s/^X//" >'PolyScan/poly_clip.c' <<'END_OF_FILE'
X/*
X * Generic Convex Polygon Scan Conversion and Clipping
X * by Paul Heckbert
X * from "Graphics Gems", Academic Press, 1990
X */
X
X/*
X * poly_clip.c: homogeneous 3-D convex polygon clipper
X *
X * Paul Heckbert    1985, Dec 1989
X */
X
X#include "poly.h"
X
X#define SWAP(a, b, temp)    {temp = a; a = b; b = temp;}
X#define COORD(vert, i) ((double *)(vert))[i]
X
X#define CLIP_AND_SWAP(elem, sign, k, p, q, r) { \
X    poly_clip_to_halfspace(p, q, &v->elem-(double *)v, sign, sign*k); \
X    if (q->n==0) {p1->n = 0; return POLY_CLIP_OUT;} \
X    SWAP(p, q, r); \
X}
X
X/*
X * poly_clip_to_box: Clip the convex polygon p1 to the screen space box
X * using the homogeneous screen coordinates (sx, sy, sz, sw) of each vertex,
X * testing if v->sx/v->sw > box->x0 and v->sx/v->sw < box->x1,
X * and similar tests for y and z, for each vertex v of the polygon.
X * If polygon is entirely inside box, then POLY_CLIP_IN is returned.
X * If polygon is entirely outside box, then POLY_CLIP_OUT is returned.
X * Otherwise, if the polygon is cut by the box, p1 is modified and
X * POLY_CLIP_PARTIAL is returned.
X */
X
Xint poly_clip_to_box(p1, box)
Xregister Poly *p1;
Xregister Poly_box *box;
X{
X    int x0out = 0, x1out = 0, y0out = 0, y1out = 0, z0out = 0, z1out = 0;
X    register int i;
X    register Poly_vert *v;
X    Poly p2, *p, *q, *r;
X
X    /* count vertices "outside" with respect to each of the six planes */
X    for (v=p1->vert, i=p1->n; i>0; i--, v++) {
X    if (v->sx < box->x0*v->sw) x0out++;    /* out on left */
X    if (v->sx > box->x1*v->sw) x1out++;    /* out on right */
X    if (v->sy < box->y0*v->sw) y0out++;    /* out on top */
X    if (v->sy > box->y1*v->sw) y1out++;    /* out on bottom */
X    if (v->sz < box->z0*v->sw) z0out++;    /* out on near */
X    if (v->sz > box->z1*v->sw) z1out++;    /* out on far */
X    }
X
X    /* check if all vertices inside */
X    if (x0out+x1out+y0out+y1out+z0out+z1out == 0) return POLY_CLIP_IN;
X
X    /* check if all vertices are "outside" any of the six planes */
X    if (x0out==p1->n || x1out==p1->n || y0out==p1->n ||
X    y1out==p1->n || z0out==p1->n || z1out==p1->n) {
X        p1->n = 0;
X        return POLY_CLIP_OUT;
X    }
X
X    /*
X     * now clip against each of the planes that might cut the polygon,
X     * at each step toggling between polygons p1 and p2
X     */
X    p = p1;
X    q = &p2;
X    if (x0out) CLIP_AND_SWAP(sx, -1., box->x0, p, q, r);
X    if (x1out) CLIP_AND_SWAP(sx,  1., box->x1, p, q, r);
X    if (y0out) CLIP_AND_SWAP(sy, -1., box->y0, p, q, r);
X    if (y1out) CLIP_AND_SWAP(sy,  1., box->y1, p, q, r);
X    if (z0out) CLIP_AND_SWAP(sz, -1., box->z0, p, q, r);
X    if (z1out) CLIP_AND_SWAP(sz,  1., box->z1, p, q, r);
X
X    /* if result ended up in p2 then copy it to p1 */
X    if (p==&p2)
X    bcopy(&p2, p1, sizeof(Poly)-(POLY_NMAX-p2.n)*sizeof(Poly_vert));
X    return POLY_CLIP_PARTIAL;
X}
X
X/*
X * poly_clip_to_halfspace: clip convex polygon p against a plane,
X * copying the portion satisfying sign*s[index] < k*sw into q,
X * where s is a Poly_vert* cast as a double*.
X * index is an index into the array of doubles at each vertex, such that
X * s[index] is sx, sy, or sz (screen space x, y, or z).
X * Thus, to clip against xmin, use
X *    poly_clip_to_halfspace(p, q, XINDEX, -1., -xmin);
X * and to clip against xmax, use
X *    poly_clip_to_halfspace(p, q, XINDEX,  1.,  xmax);
X */
X
Xvoid poly_clip_to_halfspace(p, q, index, sign, k)
XPoly *p, *q;
Xregister int index;
Xdouble sign, k;
X{
X    register int m;
X    register double *up, *vp, *wp;
X    register Poly_vert *v;
X    int i;
X    Poly_vert *u;
X    double t, tu, tv;
X
X    q->n = 0;
X    q->mask = p->mask;
X
X    /* start with u=vert[n-1], v=vert[0] */
X    u = &p->vert[p->n-1];
X    tu = sign*COORD(u, index) - u->sw*k;
X    for (v= &p->vert[0], i=p->n; i>0; i--, u=v, tu=tv, v++) {
X    /* on old polygon (p), u is previous vertex, v is current vertex */
X    /* tv is negative if vertex v is in */
X    tv = sign*COORD(v, index) - v->sw*k;
X    if (tu<=0. ^ tv<=0.) {
X        /* edge crosses plane; add intersection point to q */
X        t = tu/(tu-tv);
X        up = (double *)u;
X        vp = (double *)v;
X        wp = (double *)&q->vert[q->n];
X        for (m=p->mask; m!=0; m>>=1, up++, vp++, wp++)
X        if (m&1) *wp = *up+t*(*vp-*up);
X        q->n++;
X    }
X    if (tv<=0.)        /* vertex v is in, copy it to q */
X        q->vert[q->n++] = *v;
X    }
X}
END_OF_FILE
if test 4217 -ne `wc -c <'PolyScan/poly_clip.c'`; then
    echo shar: \"'PolyScan/poly_clip.c'\" unpacked with wrong size!
fi
# end of 'PolyScan/poly_clip.c'
fi
if test -f 'PolyScan/poly_scan.c' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'PolyScan/poly_scan.c'\"
else
echo shar: Extracting \"'PolyScan/poly_scan.c'\" \(4638 characters\)
sed "s/^X//" >'PolyScan/poly_scan.c' <<'END_OF_FILE'
X/*
X * Generic Convex Polygon Scan Conversion and Clipping
X * by Paul Heckbert
X * from "Graphics Gems", Academic Press, 1990
X */
X
X/*
X * poly_scan.c: point-sampled scan conversion of convex polygons
X *
X * Paul Heckbert    1985, Dec 1989
X */
X
X#include <stdio.h>
X#include <math.h>
X#include "poly.h"
X
X/*
X * poly_scan: Scan convert a polygon, calling pixelproc at each pixel with an
X * interpolated Poly_vert structure.  Polygon can be clockwise or ccw.
X * Polygon is clipped in 2-D to win, the screen space window.
X *
X * Scan conversion is done on the basis of Poly_vert fields sx and sy.
X * These two must always be interpolated, and only they have special meaning
X * to this code; any other fields are blindly interpolated regardless of
X * their semantics.
X *
X * The pixelproc subroutine takes the arguments:
X *
X *    pixelproc(x, y, point)
X *    int x, y;
X *    Poly_vert *point;
X *
X * All the fields of point indicated by p->mask will be valid inside pixelproc
X * except sx and sy.  If they were computed, they would have values
X * sx=x+.5 and sy=y+.5, since sampling is done at pixel centers.
X */
X
Xvoid poly_scan(p, win, pixelproc)
Xregister Poly *p;        /* polygon */
XWindow *win;            /* 2-D screen space clipping window */
Xvoid (*pixelproc)();        /* procedure called at each pixel */
X{
X    register int i, li, ri, y, ly, ry, top, rem, mask;
X    double ymin;
X    Poly_vert l, r, dl, dr;
X
X    if (p->n>POLY_NMAX) {
X    fprintf(stderr, "poly_scan: too many vertices: %d\n", p->n);
X    return;
X    }
X
X    ymin = HUGE;
X    for (i=0; i<p->n; i++)        /* find top vertex (y points down) */
X    if (p->vert[i].sy < ymin) {
X        ymin = p->vert[i].sy;
X        top = i;
X    }
X
X    li = ri = top;            /* left and right vertex indices */
X    rem = p->n;                /* number of vertices remaining */
X    y = ceil(ymin-.5);            /* current scan line */
X    ly = ry = y-1;            /* lower end of left & right edges */
X    mask = p->mask & ~POLY_MASK(sy);    /* stop interpolating screen y */
X
X    while (rem>0) {    /* scan in y, activating new edges on left & right */
X            /* as scan line passes over new vertices */
X
X    while (ly<=y && rem>0) {    /* advance left edge? */
X        rem--;
X        i = li-1;            /* step ccw down left side */
X        if (i<0) i = p->n-1;
X        incrementalize_y(&p->vert[li], &p->vert[i], &l, &dl, y, mask);
X        ly = floor(p->vert[i].sy+.5);
X        li = i;
X    }
X    while (ry<=y && rem>0) {    /* advance right edge? */
X        rem--;
X        i = ri+1;            /* step cw down right edge */
X        if (i>=p->n) i = 0;
X        incrementalize_y(&p->vert[ri], &p->vert[i], &r, &dr, y, mask);
X        ry = floor(p->vert[i].sy+.5);
X        ri = i;
X    }
X
X    while (y<ly && y<ry) {        /* do scanlines till end of l or r edge */
X        if (y>=win->y0 && y<=win->y1)
X        if (l.sx<=r.sx) scanline(y, &l, &r, win, pixelproc, mask);
X        else        scanline(y, &r, &l, win, pixelproc, mask);
X        y++;
X        increment(&l, &dl, mask);
X        increment(&r, &dr, mask);
X    }
X    }
X}
X
X/* scanline: output scanline by sampling polygon at Y=y+.5 */
X
Xstatic scanline(y, l, r, win, pixelproc, mask)
Xint y, mask;
XPoly_vert *l, *r;
XWindow *win;
Xvoid (*pixelproc)();
X{
X    int x, lx, rx;
X    Poly_vert p, dp;
X
X    mask &= ~POLY_MASK(sx);        /* stop interpolating screen x */
X    lx = ceil(l->sx-.5);
X    if (lx<win->x0) lx = win->x0;
X    rx = floor(r->sx-.5);
X    if (rx>win->x1) rx = win->x1;
X    if (lx>rx) return;
X    incrementalize_x(l, r, &p, &dp, lx, mask);
X    for (x=lx; x<=rx; x++) {        /* scan in x, generating pixels */
X    (*pixelproc)(x, y, &p);
X    increment(&p, &dp, mask);
X    }
X}
X
X/*
X * incrementalize_y: put intersection of line Y=y+.5 with edge between points
X * p1 and p2 in p, put change with respect to y in dp
X */
X
Xstatic incrementalize_y(p1, p2, p, dp, y, mask)
Xregister double *p1, *p2, *p, *dp;
Xregister int mask;
Xint y;
X{
X    double dy, frac;
X
X    dy = ((Poly_vert *)p2)->sy - ((Poly_vert *)p1)->sy;
X    if (dy==0.) dy = 1.;
X    frac = y+.5 - ((Poly_vert *)p1)->sy;
X
X    for (; mask!=0; mask>>=1, p1++, p2++, p++, dp++)
X    if (mask&1) {
X        *dp = (*p2-*p1)/dy;
X        *p = *p1+*dp*frac;
X    }
X}
X
X/*
X * incrementalize_x: put intersection of line X=x+.5 with edge between points
X * p1 and p2 in p, put change with respect to x in dp
X */
X
Xstatic incrementalize_x(p1, p2, p, dp, x, mask)
Xregister double *p1, *p2, *p, *dp;
Xregister int mask;
Xint x;
X{
X    double dx, frac;
X
X    dx = ((Poly_vert *)p2)->sx - ((Poly_vert *)p1)->sx;
X    if (dx==0.) dx = 1.;
X    frac = x+.5 - ((Poly_vert *)p1)->sx;
X
X    for (; mask!=0; mask>>=1, p1++, p2++, p++, dp++)
X    if (mask&1) {
X        *dp = (*p2-*p1)/dx;
X        *p = *p1+*dp*frac;
X    }
X}
X
Xstatic increment(p, dp, mask)
Xregister double *p, *dp;
Xregister int mask;
X{
X    for (; mask!=0; mask>>=1, p++, dp++)
X    if (mask&1)
X        *p += *dp;
X}
END_OF_FILE
if test 4638 -ne `wc -c <'PolyScan/poly_scan.c'`; then
    echo shar: \"'PolyScan/poly_scan.c'\" unpacked with wrong size!
fi
# end of 'PolyScan/poly_scan.c'
fi
if test -f 'Roots3And4.c' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'Roots3And4.c'\"
else
echo shar: Extracting \"'Roots3And4.c'\" \(4127 characters\)
sed "s/^X//" >'Roots3And4.c' <<'END_OF_FILE'
X/*
X *  Roots3And4.c
X *
X *  Utility functions to find cubic and quartic roots,
X *  coefficients are passed like this:
X *
X *      c[0] + c[1]*x + c[2]*x^2 + c[3]*x^3 + c[4]*x^4 = 0
X *
X *  The functions return the number of non-complex roots and
X *  put the values into the s array.
X *
X *  Author:         Jochen Schwarze ([e-mail])
X *
X *  Jan 26, 1990    Version for Graphics Gems
X *  Oct 11, 1990    Fixed sign problem for negative q's in SolveQuartic
X *                  (reported by Mark Podlipec),
X *                  Old-style function definitions,
X *                  IsZero() as a macro
X */
X
X#include    <math.h>
X
X/* epsilon surrounding for near zero values */
X
X#define     EQN_EPS     1e-9
X#define        IsZero(x)    ((x) > -EQN_EPS && (x) < EQN_EPS)
X
X
Xint SolveQuadric(c, s)
X    double c[ 3 ];
X    double s[ 2 ];
X{
X    double p, q, D;
X
X    /* normal form: x^2 + px + q = 0 */
X
X    p = c[ 1 ] / (2 * c[ 2 ]);
X    q = c[ 0 ] / c[ 2 ];
X
X    D = p * p - q;
X
X    if (IsZero(D))
X    {
X    s[ 0 ] = - p;
X    return 1;
X    }
X    else if (D < 0)
X    {
X    return 0;
X    }
X    else if (D > 0)
X    {
X    double sqrt_D = sqrt(D);
X
X    s[ 0 ] =   sqrt_D - p;
X    s[ 1 ] = - sqrt_D - p;
X    return 2;
X    }
X}
X
X
Xint SolveCubic(c, s)
X    double c[ 4 ];
X    double s[ 3 ];
X{
X    int     i, num;
X    double  sub;
X    double  A, B, C;
X    double  sq_A, p, q;
X    double  cb_p, D;
X
X    /* normal form: x^3 + Ax^2 + Bx + C = 0 */
X
X    A = c[ 2 ] / c[ 3 ];
X    B = c[ 1 ] / c[ 3 ];
X    C = c[ 0 ] / c[ 3 ];
X
X    /*  substitute x = y - A/3 to eliminate quadric term:
X    x^3 +px + q = 0 */
X
X    sq_A = A * A;
X    p = 1.0/3 * (- 1.0/3 * sq_A + B);
X    q = 1.0/2 * (2.0/27 * A * sq_A - 1.0/3 * A * B + C);
X
X    /* use Cardano's formula */
X
X    cb_p = p * p * p;
X    D = q * q + cb_p;
X
X    if (IsZero(D))
X    {
X    if (IsZero(q)) /* one triple solution */
X    {
X        s[ 0 ] = 0;
X        num = 1;
X    }
X    else /* one single and one double solution */
X    {
X        double u = cbrt(-q);
X        s[ 0 ] = 2 * u;
X        s[ 1 ] = - u;
X        num = 2;
X    }
X    }
X    else if (D < 0) /* Casus irreducibilis: three real solutions */
X    {
X    double phi = 1.0/3 * acos(-q / sqrt(-cb_p));
X    double t = 2 * sqrt(-p);
X
X    s[ 0 ] =   t * cos(phi);
X    s[ 1 ] = - t * cos(phi + M_PI / 3);
X    s[ 2 ] = - t * cos(phi - M_PI / 3);
X    num = 3;
X    }
X    else /* one real solution */
X    {
X    double sqrt_D = sqrt(D);
X    double u = cbrt(sqrt_D - q);
X    double v = - cbrt(sqrt_D + q);
X
X    s[ 0 ] = u + v;
X    num = 1;
X    }
X
X    /* resubstitute */
X
X    sub = 1.0/3 * A;
X
X    for (i = 0; i < num; ++i)
X    s[ i ] -= sub;
X
X    return num;
X}
X
X
Xint SolveQuartic(c, s)
X    double c[ 5 ];
X    double s[ 4 ];
X{
X    double  coeffs[ 4 ];
X    double  z, u, v, sub;
X    double  A, B, C, D;
X    double  sq_A, p, q, r;
X    int     i, num;
X
X    /* normal form: x^4 + Ax^3 + Bx^2 + Cx + D = 0 */
X
X    A = c[ 3 ] / c[ 4 ];
X    B = c[ 2 ] / c[ 4 ];
X    C = c[ 1 ] / c[ 4 ];
X    D = c[ 0 ] / c[ 4 ];
X
X    /*  substitute x = y - A/4 to eliminate cubic term:
X    x^4 + px^2 + qx + r = 0 */
X
X    sq_A = A * A;
X    p = - 3.0/8 * sq_A + B;
X    q = 1.0/8 * sq_A * A - 1.0/2 * A * B + C;
X    r = - 3.0/256*sq_A*sq_A + 1.0/16*sq_A*B - 1.0/4*A*C + D;
X
X    if (IsZero(r))
X    {
X    /* no absolute term: y(y^3 + py + q) = 0 */
X
X    coeffs[ 0 ] = q;
X    coeffs[ 1 ] = p;
X    coeffs[ 2 ] = 0;
X    coeffs[ 3 ] = 1;
X
X    num = SolveCubic(coeffs, s);
X
X    s[ num++ ] = 0;
X    }
X    else
X    {
X    /* solve the resolvent cubic ... */
X
X    coeffs[ 0 ] = 1.0/2 * r * p - 1.0/8 * q * q;
X    coeffs[ 1 ] = - r;
X    coeffs[ 2 ] = - 1.0/2 * p;
X    coeffs[ 3 ] = 1;
X
X    (void) SolveCubic(coeffs, s);
X
X    /* ... and take the one real solution ... */
X
X    z = s[ 0 ];
X
X    /* ... to build two quadric equations */
X
X    u = z * z - r;
X    v = 2 * z - p;
X
X    if (IsZero(u))
X        u = 0;
X    else if (u > 0)
X        u = sqrt(u);
X    else
X        return 0;
X
X    if (IsZero(v))
X        v = 0;
X    else if (v > 0)
X        v = sqrt(v);
X    else
X        return 0;
X
X    coeffs[ 0 ] = z - u;
X    coeffs[ 1 ] = q < 0 ? -v : v;
X    coeffs[ 2 ] = 1;
X
X    num = SolveQuadric(coeffs, s);
X
X    coeffs[ 0 ]= z + u;
X    coeffs[ 1 ] = q < 0 ? v : -v;
X    coeffs[ 2 ] = 1;
X
X    num += SolveQuadric(coeffs, s + num);
X    }
X
X    /* resubstitute */
X
X    sub = 1.0/4 * A;
X
X    for (i = 0; i < num; ++i)
X    s[ i ] -= sub;
X
X    return num;
X}
X
END_OF_FILE
if test 4127 -ne `wc -c <'Roots3And4.c'`; then
    echo shar: \"'Roots3And4.c'\" unpacked with wrong size!
fi
# end of 'Roots3And4.c'
fi
echo shar: End of archive 3 \(of 5\).
cp /dev/null ark3isdone
MISSING=""
for I in 1 2 3 4 5 ; do
    if test ! -f ark${I}isdone ; then
    MISSING="${MISSING} ${I}"
    fi
done
if test "${MISSING}" = "" ; then
    echo You have unpacked all 5 archives.
    rm -f ark[1-9]isdone
else
    echo You still need to unpack the following archives:
    echo "        " ${MISSING}
fi
##  End of shell archive.
exit 0